#!/usr/bin/perl
# 2014/11/13 頻度が極端に多かったり少なかったりするFragmentを除く

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;


if(@ARGV != 6 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [map file] -l [list of abnormal fragment] -o [output map file]\n";
}

#---------------------------------------
# map format
#---------------------------------------
# id chr1 loc1 direction1 quality1 hind1 hind_loc1 chr2 loc2 direction2 quality2 hind2 hind_loc2 type(unique or repeat) number(alignable locations) locations

my %opt;
getopts("i:l:o:", \%opt);
my $FILE_map = $opt{i};
my $FILE_list = $opt{l};
my $FILE_out = $opt{o};


#---------------------------------------
# read list file
#---------------------------------------
my %abnormal;
my $fh_list = IO::File->new($FILE_list) or die "cannot open $FILE_list: $!";
while($_ = $fh_list->getline()){
	s/\r?\n//;
	$abnormal{"$_"} = 1;
}
$fh_list->close();






my $COUNTER_passed = 0;
my $COUNTER_all = 0;

my $fh_out = IO::File->new($FILE_out, 'w') or die "cannot write $FILE_out: $!";
my $fh_map = IO::File->new($FILE_map) or die "cannot open $FILE_map: $!";
while($_ = $fh_map->getline()){
	s/\r?\n//;
	my ($id, $chr1, $loc1, $direction1, $mapQ1, $hinID1, $hinLoc1, $chr2, $loc2, $direction2, $mapQ2, $hinID2, $hinLoc2, $type, $num, $locations) = split /\t/;

	my $fragNum1 = substr($hinID1, 0, -1);
	my $fragNum2 = substr($hinID2, 0, -1);

	# 左側のHindIIIのIDとするので、右側のHindIIIに近い場合は-1する
	if($hinLoc1 > $loc1){
		$fragNum1--;
	}
	if($hinLoc2 > $loc2){
		$fragNum2--;
	}

	$COUNTER_all++;

	if((exists $abnormal{$fragNum1}) or (exists $abnormal{$fragNum2})){
		next;
	}

	$COUNTER_passed++;
	$fh_out->print("$_\n");

}
$fh_map->close();
$fh_out->close();


print "Abnormal frequency fragment filter\n";
printf "passed reads: %d (%.3f)\n", $COUNTER_passed, $COUNTER_passed / $COUNTER_all * 100;
printf "All readds: %d\n", $COUNTER_all;
